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Abstract: We generalize and test the recent "ab initio" self-consistent 
(AISC) time-independent semiclassical laser theory. This self-consistent 
formalism generates all the stationary lasing properties in the multimode 
regime (frequencies, thresholds, internal and external fields, output power 
and emission pattern) from simple inputs: the dielectric function of the 
passive cavity, the atomic transition frequency, and the transverse relaxation 
time of the lasing transition. We find that the theory gives excellent quantita- 
tive agreement with full time-dependent simulations of the Maxwell-Bloch 
equations after it has been generalized to drop the slowly-varying envelope 
approximation. The theory is infinite order in the non-linear hole-burning 
interaction; the widely used third order approximation is shown to fail 
badly. 
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The Maxwell-Bloch (MB) equations provide the foundation of semiclassical laser theory U 
and are the simplest description which captures the full space-dependent non-linear behavior 
of a laser. These time-dependent equations can be simulated to determine the stationary lasing 



state. However time-independent methods to find these stationary properties in the multi-mode 
regime for an open laser cavity did not exist until the recent work of Tureci et al. J21 [3] 0] 
presented an "ab initio" self-consistent (AISC) formalism which generates all of the lasing 
properties including the output power and emission pattern from a few simple inputs. The 
laser cavity can be of arbitrary complexity and openness, including, e.g., chaotic dielectric disk 
lasers Q, photonic lattice defect mode lasers Q and random lasers fl4j[7). Here we generalize 
this infinite order non-linear theory by extending it beyond the slowly-varying envelope ap- 
proximation. With this improvement it gives remarkably good agreement with time-dependent 
simulations of the Maxwell-Bloch (MB) equations, while the standard third order approxima- 
tion to the non-linear hole-burning interaction fails badly. 

The AISC theory builds on the original ideas of Haken and coworkers [H[9) that the inversion 
of the lasing medium will be approximately time-independent when <C y± , A (where y± is the 
transverse (polarization) relaxation rate, A is the typical mode spacing, and y\\ is the longitudinal 
(inversion) relaxation rate). The only significant approximation in the theory presented below 
is this approximation of stationary inversion (SIA) (well-satisfied in many lasers of interest). In 
addition to the excellent agreement we find between the AISC theory and the MB simulations 
when the ratios yji/y_L,yji/A are very small, we develop below a perturbative treatment of the 
beating terms which are neglected in SIA to extend the theory to larger values of these ratios. 
The key improvements contained in the AISC theory are: 1) Treatment of the openness of the 
cavity exactly. 2) Inclusion of the space-dependent non-linear modal interactions (spatial hole- 
burning) to all orders, in contrast to standard third-order treatments [Q. We show below that 
the third order treatment fails quantitatively and qualitatively even for the simple laser resonator 
we study here. 

To perform a well-controlled comparison of MB and AISC results we chose to study the 
simple one-dimensional microcavity edge emitting laser 121 consisting of a perfect mirror at 
the origin and a dielectric region of uniform index «o and length L terminating abruptly on air 
(see inset, Fig. [T). The MB equations are simulated in time and space using a FDTD approach 
for the Maxwell equations, while the Bloch equations are discretized using a Crank-Nicholson 
scheme. To avoid solving a nonlinear system of equations at each spatial location and time 
step, we adopt the method proposed by Bidegaray ifTol . in which the polarization and inversion 
are spatially aligned with the electric field, but are computed at staggered times, along with the 
magnetic field. Modal intensities are computed by a Fourier transform of the electric field at 
the boundary after the simulation has reached the steady state. 

The AISC theory consists of a set of coupled non-linear time-independent integral equations 
of size equal to the number of lasing modes at a given pump. The AISC theory presented in 
Refs. [2, 3] is a solution to the MB equations [ l] after two standard simplifications, the rotating- 
wave approximation (RWA) and slowly-varying envelope approximation (SVEA). The SVEA 
involves factoring out the rapid time-dependence of the electric field and the atomic polarization 
field at the atomic frequency, k a , (here and below we set c = 1 and use frequency and wave vector 
interchangeably), and neglecting the second time derivatives in the Maxwell wave equation 
of the remaining envelope function of the fields. The resulting non-linear system contains 
only first time-derivatives of the fields and the inversion and is sometimes referred to as the 
Schrodinger-Bloch (SB) equations fTTI . The SVEA works well when the cavity frequencies 
are negligibly shifted from the atomic frequency. For microcavities these shifts need not be 
negligible; our simulated cavities have nok a L = 60 corresponding to roughly ten wavelengths 
of radiation inside the cavity, approaching the microcavity limit. For the case studied here, we 
find noticeable discrepancies between MB solutions and the AISC theory of Refs. [2-4] which 
incorporates the SVEA (see Fig. [T). This motivated us to generalize the AISC theory to drop 
the SVEA (the RWA is found to be well-satisfied in all cases). 
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Figure 1: Modal intensities as functions of the pump strength Dq in a one-dimensional mi- 
crocavity edge emitting laser of y± = 4.0 and y\\ = 0.001. (a) n = 1.5, k a L = 40. (b) n = 3, 
k a L = 20. Square data points are the result of MB time-dependent simulations; solid lines are 
the result of time-independent ab initio calculations (AI) of Eq. (Q]). Excellent agreement is 
found with no fitting parameter. Colored lines represent individual modal output intensities; 
the black lines the total output intensity. Dashed lines are results of AISC calculations when 
the slowly-varying envelope approximation is made as in Ref. [3] showing significant quantita- 
tive discrepancies. For example, in the n = 3 case the differences of the third/fourth thresholds 
between the MB and AISC approaches are 46% / 63%, respectively, but are reduced to 3% and 
15% once the SVEA is removed. The spectra at Do = 10 and the gain curve are shown as insets 
in (a) and (b) with the solid lines representing the predictions of the AISC approach (Eq. (HJ) 
and with the diamonds illustrating the height and frequency of each lasing peak. The schematic 
in (a) shows a uniform dielectric cavity with a perfect mirror on the left and a dielectric-air 
interface on the right. 



The generalization was as follows. Again stationary periodic solutions are assumed for 
the electric field E(x,t) = 2Re [e(jr,f)] = 2Re [J^ ^ (jt) exp (— ik^t)\ and for the polarization 
fields, which oscillate at unknown lasing frequencies, k^. The spatial variation of the field am- 
plitude *Pm (x) is also unknown, and not assumed to be a cavity resonance, but is determined 
self-consistently. For a high finesse cavity and near the first threshold it was shown 13 that 
(x) is well-approximated by a single cavity resonance, but above threshold and for lower 
finesse this is not at all the case OH]. The treatment of the matter equations does not involve 
the SVEA and is exactly the same as in Ref. [2] , where the key assumption is stationary in- 
version, which allows the non-linear polarization in the Maxwell equation to be replaced by 
a non-linear function of the electric field itself. The new element is that we keep the second 
time-derivative of the polarization in the Maxwell equation and evaluate it by differentiating 
the polarization equation. The resulting improved AI/MB equations for the mode functions ^ 
and the frequencies are: 

™ (r) in k l [ Do(x)G(x,x>;k,)V,(x>) 

mW Yx-Kh-K) kl J dj: e(jO(l+E v r v pF v (jO|2) ' 

Here G{x ) x';k fl ) is the Green function of the open cavity, Y v = r(k v ) is the gain profile 
evaluated at k v , Dq(x) = Dq(1 + do(x)) is the spatial pump, and e(x) = n 2 (x) is the dielec- 
tric function of the cavity (for the microcavity edge emitting laser n(x) = «o and the pump 
is assumed uniform (do(x) = 0)). Electric field and pump strength are measured in units 
e c = hy/Y±_Y\\/2g,Do c = \%k 2 a g 2 /hyx, where g is the dipole matrix element of the gain medium. 
With a slight change in notation this equation differs from that derived in Ref. [2] only by the 
additional factor k^/k 2 multiplying the integral. This is consistent with the expectation that 
the SVEA is good when the lasing frequency is very close to the atomic frequency; incorpo- 
rating this change into the iterative algorithm for solving the system is trivial, but crucial for 
quantitative agreement with the current MB simulations. However we do not find qualitative 
changes from dropping the SVEA, either for this simple laser or for the more complicated 
two-dimensional random laser studied elsewhere J4). 

In Figs, la, lb we show results for no = 1.5,3, finding remarkable agreement between MB 
and AI approaches with no fitting parameters for the case of three mode and four mode lasing 
respectively. Both thresholds, modal intensities and frequencies are found correctly by the 
AISC approach. Note that all of these quantities are direct outputs of the AISC theory, whereas 
they must be found by numerical fourier analysis of the MB outputs, which can introduce some 
numerical error. If the earlier AI approach is used with the SVEA then significant discrepancies 
are found, for example for «o = 3 the threshold of the third mode is found to be higher than 
from MB while the fourth threshold is too low. Note however that the theory with the SVEA 
does get the right number of modes and the correct linear behavior for large pumps. We believe 
that this original AISC approach does solve very accurately the Schrodinger-Bloch equations 
and that the same discrepancies would arise between MB and SB simulations, although we 
haven't confirmed this. 

Almost all studies of the MB equations in the multimode regime have used the approximation 
of treating the non-linear modal interactions to third-order (near threshold approximation) and 
in fact this approximation is used quite generally and uncritically throughout laser theory. From 
examination of the form of Eq. (Q]) it is clear that it treats the non-linear interactions to all 
orders, while the third order treatment would arise from expanding the denominator to the 
leading order in li'vl 2 - This third-order version of the AISC theory then becomes similar to 
standard treatments of Haken [8 9], with the improvement of correctly treating the openness of 
the cavity and the self-consistency of the lasing modes in space [2 3 1. An early version of this 
improved third order theory was found to have major deficiencies: it predicted too many lasing 



modes and the intensities did not scale linearly at large pump, but exhibited a spurious saturation 
(21 [l2l. In Fig. [2]we present comparisons of the third order approximation to Eq. (HJ, which is 
improved over Ref. [2] because it drops the SVEA and the "single-pole approximation" used 
there. We find that this improved third order theory still does a very poor job of reproducing the 
multimode MB results: it still predicts too many modes (in this case seven, when there should 
only be four at Do = 10 in Fig. [2]), and shows the same spurious saturation as found earlier J2) 
because the third-order approximation cannot give the correct linear behavior for large pumps. 
The infinite order treatment of Eq. fl} is both qualitatively and quantitatively essential. 




Pump strength D Q / D Qc 

Figure 2: Modal intensities as functions of the pump strength Do in a one-dimensional micro- 
cavity edge emitting laser of n = 3, k u L = 20, y± = 4.0 and y| = 0.001; the solid lines and data 
points are the same as in Fig. lb. The dashed lines are the results of the third order approx- 
imation to Eq. ((TJ. The frequently used third order approximation is seen to fail badly at a 
pump level roughly twice the first threshold value, exhibiting a spurious saturation not present 
in the actual MB solutions or the AI theory. In addition, the third order approximation predicts 
too many lasing modes at larger pump strength. For example, it predicts seven lasing modes at 
Do = 10, while both the MB and AISC show only four. Right inset just shows the same data on 
a larger vertical scale. 

The central approximation required for Eq. ((TJ is that of stationary inversion. Previous 
work by Haken J9| argued that SIA holds for the MB equations when y| <C y±,A , where 
Y\\ is the relaxation rate of the inversion and A is the frequency difference of lasing modes. 
For a typical semiconductor laser y ± — 10 12 — 10 13 s~' and y| ~ 10 8 — I0 9 s~ l , or equivalently 
)j| /yj_ ~ 10~ 3 — 10~ 5 . For the microcavity edge-emitting laser we are modeling we took y±_ , A ~ 
1 and y« = 10~ 3 , and we found the excellent agreement shown in Figs, la, lb. In addition, direct 
analysis of the inversion vs. time obtained from the MB simulations confirm very weak time- 
dependence in the steady state, justifying the use of SIA. The previous work did not develop 
a systematic theory in which yj|/yj_,yj|/A appear as small parameters in the lasing equations, 
allowing perturbative treatment of corrections to the SIA; we are now able to do this within the 
AI formalism. 

Note first that in Eq. ([TJ the electric field is measured in units e c ~ ^/yybut, unlike yx, Y\\ does 
not appear explicitly. Hence the solutions of Eq. (fTJ) depend on yji only through this scale factor, 
and Eq. ([TJ makes the strong prediction of a universal overall scaling of the field intensities: 
\E(x) \ 2 ~ Y\\ when dimensions are restored. The perturbative corrections to Eq. (UJ are obtained 



by including the leading effects of the beating terms between the different lasing modes which 
lead to time-dependence of the inversion at multiples of the beat frequencies. These population 
oscillations non-linearly mix with the electric field and polarization to generate all harmonics 
of the beat frequencies in principle, but the multimode approximation assumes all the newly 
generated fourier components of the fields are negligible. The leading correction to this ap- 
proximation is to evaluate the effect of the lowest sidebands of population oscillation on the 
polarization at the lasing frequencies and on the static part of the inversion, both of which will 
enter Eq. (HJ. For simplicity we present a sketch of the correction to Eq. ([T) in the two-mode 
regime; details and the straightforward generalization to more modes will be given elsewhere. 
We start by writing the electric field and the polarization in the multiperiodic forms e(x,t) = 



2Re 



L^=i l IV(*)exp(-;V) and p{x,t) = 2Re l^ =1 p^(x)exp(-ik^t 



-2 



and allow for the 



first two side-bands at the beat frequency A = k\ — k 2 in the inversion, so that the total inver- 
sion is D(x,t) =D s (x)+d+(x)sxp(—iAt)+d-(x)exp(+iAt) where the real quantity D s (x) is 
the time-independent part of the inversion and d+(x) = d~(x)*. This ansatz is inserted into 
the Bloch equation for the inversion 0; solving for the component of the inversion equation 
which oscillates at exp (—/A?) relates d+ (x) to the product of the field and polarization and then 
substitution of the zeroth order result for the polarization in terms of the zeroth order static 
inversion D s gives 

d+[x) = I [Vl (fl-yf l] = jmM)Df\x)^(x)^* 2 (x)/e 2 c (2) 

where the dimensionless function f(k\,k 2 ) = — (i + A/ (2y±))/ (1 + k\k 2 — iA/y±),k v = (k v — 
k a ) /y_j_ and the fields are not yet measured in units of e c . The component d + will mix with the 
field ^2 to yield a contribution to the polarization p\ at frequency k\, ( and similarly aL and 
*Fi mix to contribute to p 2 ), 

(1) ll + ^/K)f(k lM )\V 2 {x)\ 2 /el] (0) 
P\ V) = ¥l y± - iikl -k a ) °> (3) 

where p 2 \x) is obtained by interchanging subscripts 1 , 2. The correction to the AI formalism 
is the term in the numerator explicitly proportional to the small parameter 7[|/A. However, 
having found a correction to the polarizations pi,p2 we must then evaluate their contribution 
to the static inversion. We find 

D ( X ) = £° ( 4 ) 

lW l + L v r v |^ v W|2 + (y | /y ± )^ l! fc 2 )|^ 1 (^ 2 (x)|2' (> 

where the dimensionless function g(k\ 7 k 2 ) = (2 + kf +^|)(1 — kik,2)/[(l +k\){\ +^2)] 2 an d 
now and below electric fields have been scaled by e c . Note that the correction term in Eq. (01 is 
explicitly proportional to the second small parameter, yi\/y±. For our simulations ji_ ~ A and 
the functions f(k\,hi),g{ki,k%) are order unity. The full correction to the non-linear polariza- 
tion in Eq. (01 is obtained by replacing Dj with D s of Eq. ©. The corrected polarization 
leads to corrected version of Eq. ([]]) of the AISC theory: 

iyJh k\ I dx> (l + r {f(k l ,k 2 )\V l (x')\ 2 )G(x,x';k 2 )V2(x') 
n-i(k 2 -k a )k 2 J e(xf) (l+I v r v |1' v (x')| 2 + ^(^i^2)| v I , i(^) l I , 2(^)| 2 )' 



(x) satisfies the same equation with the subscripts 1 and 2 interchanged. 



Eq. © predicts corrections to the universal behavior, ^(jc)! 2 ~ yh, found in Fig. Q] There 
is no correction to the first mode below the second threshold as the correction terms all vanish 
(there needs to be two modes to have beats). However, the theory predicts a non-trivial correc- 
tion to the threshold of the second mode. Note that the correction to the numerator in Eq. (O 
does not vanish below the second threshold but contributes self-consistently to that threshold. 
This correction can be regarded as modifying the dielectric function of the microcavity to take 
the form e'(x) = e(x)/[l + ^f(k2,ki)\ x Vi(x)\ 2 ]; the effective dielectric function then becomes 
complex and varying in space according to the intensity of the first mode. This in turn changes 
the threshold for the second mode. If the modes are on opposite sides of the atomic frequency 
and ki < k\, the imaginary part of the effective index is always amplifying and tends to de- 
crease the second threshold; we find this effect dominates over the change in the real part and 
increasing y\\ uniformly decreases the thresholds. The opposite effect is possible and observed 
in other cases we have studied (not shown). In Fig. [3]we show the results of MB simulations as 
7j| is varied from 0.001 to 0. 1 (with y± = 4). Note that the universal behavior (in units scaled by 
is well obeyed until Y\\ = 0.1, encompassing most lasers of interest. The qualitative effect 
predicted by the perturbation theory is clearly seen, the higher thresholds are reduced as y\\ is 
increased. The effect is small for the 2nd threshold but large for the third as we expect as the 
corrections scale with the product of the intensities of lower modes. The inset to Fig. [3]shows 
that the perturbation theory for the third threshold (a suitable generalization of the two-mode 
Eq. (f5]l) yields semiquantitative agreement with the threshold shifts found from the simulations. 
Detailed comparisons between multimode perturbation theory and simulations above threshold 
will be given elsewhere. Note that the simulations also find additional modes turning on for 
J\\ =0.1 but their intensities are very small and they are not shown in Fig. [3] 
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Figure 3: Modal intensities for the microcavity edge emitting laser of Fig. la as y\\ is varied 
(n = 1.5, k a L = 20, /x = 4.0 and mode-spacing Aw 1.8. Solid lines are AISC results from 
Fig. la (with yj| = 0.001); dashed lines are y\\ = 0.01 and dot-dashed are y\\ = 0.1. The color 
scheme is the same as in Fig. la. The inset shows the shifts of the third threshold as a function of 
Y\\ ■ The perturbation theory (squares, with the line to guide the eyes) predicts semiquantitatively 
the decrease of the threshold as y» increases found in the MB simulations. The MB threshold 
is not sharp and we add an error bar to denote the size of the transition region. 

In conclusion, for the case studied, the recently developed ab initio self-consistent laser the- 
ory without the slowly varying envelope approximation presents a very accurate solution of the 
steady-state Maxwell-Bloch semiclassical lasing equations without solving the time-dependent 



problem. Third order treatments fail badly and our infinite order treatment is essential. The 
theory is a well-controlled expansion in the small parameters 7||/A, Y\\/Y± an d leading correc- 
tions in these small parameters can be evaluated and understood qualitatively. The AI method 
is flexible and more convenient than MB simulations for complex two-dimensional and three- 
dimensional microcavities and provides more physical insight into the lasing properties. The 
accuracy of the method suggests it can be useful in the analysis and design of novel laser sys- 
tems. 
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